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Abstract - In this paper we give an introduction to the Boltzmann equation for neutrino transport 
used in core collapse supernova models as well as a detailed mathematical description of the Isotropic Diffu- 
sion Source Approximation (IDSA) established in [6]. Furthermore, we present a numerical treatment of a 
reduced Boltzmann model problem based on time splitting and finite volumes and revise the discretization 
of the IDSA in |6] for this problem. Discretization error studies carried out on the reduced Boltzmann model 
problem and on the IDSA show that the errors are of order one in both cases. By a numerical example, a 
detailed comparison of the reduced model and the IDSA is carried out and interpreted. For this example the 
IDSA modeling error with respect to the reduced Boltzmann model is numerically determined and localized. 

Resume - Dans cet article, nous donnons une introduction a Pequation dc Boltzmann pour le trans- 
port des neutrinos dans les modeles de supernovae a effondrement de coeur ainsi qu'une description detaillee 
de 1' Isotropic Diffusion Source Approximation (IDSA) developpee dans [5]. De plus, nous presentons le traite- 
ment numerique d'un modele de Boltzmann simplifie base sur une decomposition en temps de l'operateur et 
sur un algorithme de volumes finis ainsi que l'adaptation de la discretisation de l'IDSA de [6| a notre modele. 
Les etudes de l'erreur de discretisation faites sur le modele de Boltzmann simplifie et sur l'IDSA montrent 
que les erreurs sont d'ordre un dans les deux cas. A l'aide d'un exemple numerique, nous comparons et 
interpretons en detail les deux modeles. Pour cet exemple, l'erreur de modelisation de l'IDSA par rapport 
au modele de Boltzmann simplifie est determinee numeriquement et localisee. 

Introduction 

Modelling neutrino transport is crucial for the simulation of core collapse supernovae since more 
than 99% of the released gravitational binding energy is carried away by neutrinos [HI p. 361] that 
are also assumed to feed the shock leading to the explosion. However, full 3D Boltzmann neutrino 
transport models are still computationally too costly to solve. The Isotropic Diffusion Source 
Approximation (IDSA) intends to capture the most important processes of neutrino transport 
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while being computationally feasible [6j. The main idea of the IDSA is to consider an additive 
decomposition of the neutrino distribution function into a trapped and a free streaming particle 
component. The resulting equations related to these two components are supposed to be coupled 
by a source term. Both the source term and the equations which are reduced to the main physical 
properties of the two particle components are derived. The source term is based on a diffusion limit 
and, for non-diffusive regimes, limited from above and below on the basis of the free streaming and 
reaction limits. 

In the first part of this paper, we give an introduction to the 0(v/c) Boltzmann equation for 
radiative transfer in the comoving frame (Section[T]) as well as to the IDSA of this model (Section[2]). 
This part is based on the presentation in [6], however, it is more mathematically oriented and, in 
particular, the derivation of the diffusion source, which is only sketched in [6], is presented in a 
comprehensive manner. As in [6], we restrict ourselves to the spherically symmetric case here, 
however, both models can be extended to the nonsymmetric 3D case. 

In the second part of the paper, we introduce a reduced Boltzmann model problem in which 
we mainly assume frozen background matter and give the corresponding IDSA (Section [3]). For 
this reduced model we present a numerical solution technique based on time splitting between the 
reaction and the transport part. While the reaction part has an analytical solution, a conservative 
formulation can be found for the transport part for which a finite volume scheme is established 
(Section [4|. Here we also recall the discretization of the IDSA given in [6] and adapt it to the 
reduced model. On this basis, we present several numerical results (Section [5]). To start with, 
convergence studies are carried out on the reduced Boltzmann model problem and on the IDSA in 
order to have realistic estimates of the discretization error. The results show that the errors are of 
order one in both cases. Then we establish a numerical example in order to compare the behaviour 
of solutions of the reduced model and of the IDSA. Extensive interpretations of the solutions in 
reaction and free streaming regimes and in the transition regime are given. For this example we also 
find a considerable IDSA modeling error with respect to the reduced Boltzmann model numerically 
and show that the error is mainly localized in the transition regime. 



1 Boltzmann's radiative transfer equation 

The 0(v/c) Boltzmann's equation is a widely accepted model for the radiative transfer of neutrinos 
in core collapse supernovae, see, e.g., [7J El E] . We present the version used in [B] and write it in 
short as 

ldf df „ df „ df 
c at or Ofi ouj 



+ ^ + F ^. + F "7T.= J-Xf + C(f) . (1.1) 



— " JU) 

T>V) 



Equation (1.1) represents the special relativistic transport equation for massless fermions up to 



the order 0(v/c), i.e., the neutrinos are considered to move with light speed c while the background 



matter moves with velocity v. We refer to [8, §95] for a derivation of (1.1) and [9] for a discussion 
of possibly additional 0(v/c) terms that might have to be considered. 
The equation describes the evolution of the distribution function 

/: [0,7] x (0,12] x [-1,1] x (0,£]-> [0,1], (t,r,^w) H> f(t,r,fj,,u) , (1.2) 

of the neutrinos. Here, T > is some end time, R > some maximal radius and E > some 



maximal neutrino energy. We consider ( 1.1 ) to be given in spherical symmetry, i.e., (0, R] represents 
the spatial domain Q, C M 3 , which is the ball around the origin with radius R. The variable 
H G [— 1, 1] is the cosine of the angle between the outward radial direction and the direction of 
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neutrino propagation and uj E R + is the neutrino energy, both (the phase space variables) given in 
the frame comoving with the background matter, cf. p. 640]. 

With regard to the notation we use the Lagrangian time derivative 

a[ _ d j9 
dt dt dr ' 

in the comoving frame. Furthermore, we have 

where F®j^- accounts for the change in propagation direction due to inward or outward movement 

of the neutrino. The second term F^^- represents the aberration (i.e., the change observed in the 
comoving frame) of the neutrino propagation due to the motion of the background matter with the 
density p. Finally, the product of the factor 



9 / d In p 3v 
V cdt cr I cr 



uj . 



with ^£ accounts for the (Doppler-) shift in neutrino energy due to the motion of the matter. The 
left hand side of the Boltzmann equation is abbreviated by T>(f) where T> is a linear operator. 

The dependency of T> on the background matter occurring in the comoving frame vanishes in 
case of a static background with frozen matter where one can pass to the laboratory frame by 
setting = = 0, i.e., 

1 df df 1 2 ,9/ . _ 



and using Lorentz transformed quantities in (1.3), see \B, §90, §95] 



Although in the infall phase [2j p. 787] it is enough to consider only electron neutrinos v e , for 
postbounce simulations [fH p. 1179] one needs at least two Boltzmann equations in order to obtain 
the transport of both electron neutrinos v e and electron antineutrinos v e . In general, one needs to 
include muon and tau neutrinos and their antiparticles, too [Hj. All different types of neutrinos are 
transported independently, so that in general one has to deal with up to six Boltzmann equations 
that are, however, coupled via their right hand sides. Since all these equations have the same basic 
structure, it is enough for our purpose to consider only one Boltzmann equation as a prototype. 



The source and sink terms on the right hand side of (1.1) account for interaction of neutrinos 



with background matter. They include emission and absorption 

+ p f± n + v e , 
e + + n f± p + v e , 

of electron neutrinos v e or electron antineutrinos v e by protons p or neutrons n, the forward reactions 
being known as electron e~ or positron e + capture. Analogous reactions occur in case of electron 
e~ or positron e + capture of nuclei. They depend on the state of the background matter and result 
in neutrino emissivity j(uj) and absorptivity x( w ) whose sum is the neutrino opacity x = J + X m 



(1.1), cf. [7j and [6[ p. 1177]. Concrete formulas for j(uj) and x( w )j which are nonlinear in uj, have 



been derived in pp. 822-826] for both electron neutrino and its antiparticle. 
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By the term C(f), the right hand side of (1.1) also accounts for isoenergetic scattering of 
neutrinos (or antineutrinos) on protons, neutrons and nuclei. C(f) is a linear integral operator in 
/, the so called collision integral, and is given by 



C(f(t,r,n,u})) 



to 



c(hc) 



R(H, //, uj) (f(t, r, //, uj) - f(t, r, fi, uj)) d/jf 



[1.4) 



where the isoenergetic scattering kernel R([i, fjf, uj) is symmetric in \x and // and depends nonlinearly 
on all its entries, see pp. 806/7, 826-828] for concrete formulas that also exhibit the dependency 



of R on the background matter. In (1.4) Planck's constant is denoted by h, the term corresponding 
to f(t,r,fi',uj) accounts for in-scattering while the term corresponding to f(t,r,fj,,uj) represents 
out-scattering {6^ p. 698]. Note that if / does not depend on \x we have C(f) = 0. As an immediate 
consequence of the symmetry of the collision kernel with respect to \x and \j! we obtain for any / 



C(f) dfx = 0. 



[1.5) 



Further possible source terms stemming from neutrino interactions with the background matter 
such as, e.g., neutrino electron scattering (cf. p. 774]), are neglected in [Bj. We abbreviate the 



right hand side of Boltzmann's equation (1.1) by j + J(f) where J(f) is linear in /. 



2 Isotropic Diffusion Source Approximation (IDSA) 

In this section we give an introduction to the Isotropic Diffusion Source Approximation (IDSA) 



that has been developed in (6). The aim of this approximation of Boltzmann's equation (1.1) is to 



reduce the computational cost for its solution, making use of the fact that ( 1.1 ) is mainly governed 



by diffusion of neutrinos in the inner core and by transport of free streaming neutrinos in the 
outer layers of a star. The following ansatz for the IDSA intends to avoid the solution of the full 
Boltzmann equation in a third domain in between these two regimes as well as the detection of the 
corresponding domain boundaries. 

2.1 Ansatz: Decomposition into trapped and streaming neutrinos 

We assume a decomposition of 

f = f + r 

on the whole domain into distribution functions / and f s supposed to account for trapped and for 
streaming neutrinos, respectively. 

With this assumption and by linearity V(f) = V{f) + V(f s ) and J(f) = J(f) + J{f s ), 



solving the Boltzmann equation (1.1), i.e., 

V{f)=j + J{f), 

is equivalent to solving the two equations 

T>(f)=j + J(f)-E, 
V(f) = J(f s ) + X, 



(2.1) 

(2.2) 
(2.3) 



with an arbitrary coupling term S = T,(t,r, fj,,u, /* , f s , j, J). For the IDSA one establishes ap- 
proximations of the angular mean of these two equations arising from physical properties of 
trapped and streaming particle, respectively, and one determines an appropriate coupling func- 
tion T 1 {t,r,^uj,f\f s ,j,J). 
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2.2 Hypotheses and their consequences for trapped and streaming particle 
equations 

Concretely, one uses the following hypotheses. First, the trapped particle component /* as well as 



S are assumed to be isotropic, i.e., independent of [i. Taking the angular mean of equation (2.2) 
this immediately leads to the trapped particle equation 



df 

cdt ' 3 cdt 
in which we slightly abuse the notation 



ldlnp df 

duo 



(2.4) 



(2.5) 



concerning the domain of definition of the isotropic / given in (1.2). The isotropic source term £ 



is treated in the same way. Here, and in what follows, we always assume that we can interchange 
the integral and the differentiation operators. 

Next, /* is assumed to be in the diffusion limit, which is physically at least justified for the 
inner core of the star. In order to derive the diffusion limit, a Legendre expansion of the scattering 
kernel R(/i,/i',La) with respect to its angular dependence, truncated after the second term, is used 
in [6j App. A] for an approximation of the collision integral, see Subsection 2.3 In fact, this 



approximation is essential for the derivation of the diffusion limit and thus the corresponding 



definition of £ that we will provide in Subsection 2.4 



Second, the streaming particle component f s is assumed to be in the free streaming limit. This 



justifies to neglect the collision integral in (2.3), which by (1.5), however, vanishes anyway after 



angular integration. Furthermore, it justifies to neglect the dynamics of background matter so that 



one can use the laboratory frame formulation ( 1.3 ) of ( 1.1 ) with frozen matter (here, we also neglect 
the Lorentz transformation). For the same reason one can assume the free streaming particle to 



be in the stationary state limit and drop the time derivative in (1.3) which then, after angular 
integration, becomes the streaming particle equation 



Id 



f'/idn 



f s dfi + S . 



(2.6) 



Since in spherical symmetry, with the radial unit vector e r , the gradient and the divergence oper- 
ators are given by 



d 1 d 

W>(r) = — V(r) e r and V • F(r) = {r 2 F r (r)) 



(2.7) 



for scalar fields ip : R -> R and vector fields F = (F r ,0,0) : R -> R 3 (see, e.g., P p. 679]), this 
equation can be reformulated as a Poisson equation for a spatial potential i/j of the first angular 
moment of f s . 



For the solution of (2.6), the approximate relationship 



f s {u)ixdn 



1 + 



max(r, R v {oj)) 



f s (u)dn, 



(2.8) 



between the particle flux and the particle density of streaming neutrinos, which has been suggested 
by Bruenn in [5], is applied. Here, R v (u) > is the energy dependent radius of the neutrino scat- 
tering spheres. In addition to spherical symmetry, this approximation is based on the assumption 
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that all free streaming particles of a given energy u are emitted isotropically at their correspond- 
ing scattering sphere [6] p. 1178]. As a consequence, the flux can be expressed as a product of 



the particle density and the geometrical factor in (2.8). Note that this factor is equal to 1 when 



T < Rvi^), which follows from the isotropy of / inside the scattering spheres, and increases up 
to 2 in the limit r — > oo. The latter expresses the fact that the neutrinos tend to stream radially 
outwards so that the distribution function / accumulates at ll = 1. 



2.3 Legendre expansion of the scattering kernel 

As mentioned in the last subsection, we now seek for an approximation of the collision integral 
by a Legendre expansion of the scattering kernel. For an introduction to Legendre expansions 
by spherical harmonics we refer to j!0| pp. 302, 391-395]. Concretely, the Legendre series for 
^^R(ll,li',u)) reads 



or 



c{hc) 



r-R(/t, //, w) 



— 23(2Z + l)0,(w) / Pi(cos0)dcp, 

™ 1=0 ^° 



(2.9) 



with the Legendre polynomials Pi, I = 0, 1, . . ., where 6 is the angle between the incoming and the 
outgoing particle and 

cos(6>) = mm' + [(1 - li 2 )(1 - ll' 2 )]^ cos (p. (2.10) 
With the first two Legendre polynomials given by 



P (cos(9) = 1 , Pi (cos ( 
truncation of the series after the second term provides 



cos 9 . 



c(hc) 



T R(n,n',u) 



1 

4"7T 



2vr f-2n 

0o(w) j 1 dip + 30i (cj) / cos(9)dip 
o Jo 



Inserting this in the collision integral (1.4), without explicitly mentioning the dependency on t, r 
and cj, one obtains 



(0o + 30 ly u//) (/(//) - f(n))dn' = -0 O / + 0o 



fd/u, + 3/x0i 



1 



fiid/j,. 



(2.11) 



which is a affine function in [i expressed in terms of /, the zeroth and first angular moments of 
/ and the opacities 0o and 0i. Together with the emissivity and absorptivity of neutrinos by the 



matter on the right hand side of (1.1) the latter gives rise to the definition of the neutrino mean 
free path 

X -=~ \ - = - T -• ( 2 - 12 ) 

J + X + 00 - 01 X + 00 - 01 

This definition is motivated by the fact that A/3 occurs as the diffusion parameter in the diffusion 



limit of the Boltzmann equation that will be derived in the following subsection, see (2.25). It is 



clear that the smaller the diffusion parameter is, the smaller the diffusion of neutrinos represented 



by the diffusion term in (2.25) becomes which physically corresponds to a smaller mean free path. 



Finally, we remark that a collision kernel which can be expanded as in (2.9) with (2.10) is always 



symmetric in ll and // since cos(#) in (2.10) has this property. For the same reason, (2.11) as well 



as any truncation of the Legendre series in (2.9) is symmetric in ll and ll' 
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2.4 Derivation of the diffusion limit 



Now we outline the line of thought for the derivation of the diffusion limit given in [6] . It is based 
on the truncation of the Legendre expansion of the collision kernel after the second term given in 
the previous subsection. With (2.11) we obtain the equation 

1 



2>(/)=j-(x + 4>o)/ + <£o; 



-l z 



ffidn 



(2.13) 



-i 



as an approximation of the Boltzmann equation (1.1). The basic idea for the derivation of the 
diffusion limit is to exploit the special structure of the right hand side in (2.13), i.e., the fact that 
/ can be expressed in terms of T>(f) and the zeroth and first angular moment of /. By taking the 
zeroth and first angular moments of ( |2.13 ), one can therefore express these moments of / in terms 
of the zeroth and first angular moments of T>(f), i.e., eliminate them in (2.13) and thus express / 
solely in terms of T>{f) and its zeroth and first angular moment. 

Concretely, since j is isotropic and the collision kernel in (2.11) that appears in the right hand 
side of (2.13) is symmetric in \i and fx', i.e., (1.5) is applicable, the zeroth moment of (2.13) reduces 
to 



V{f)dn = j-x 



-i 



fdfi. 



-i 



i.e., 



fdfi 



X 



(2.14) 



(2.15) 



In the first moment of (2.13), the first and third summand on the right hand side vanish 
since they are independent of fi and the first angular moment of a constant is zero. Considering 
\ J —1 3[i 2 dfi = 1, the coefficient in the last summand reduces to <j)\ so that we obtain 



V{f)fidfi = -(x + <t>o~ <t>i): 



-i 



ffidfi, 



-l 



which, by (2.12), leads to 



ffidfi 



(2.16) 



-i 



With (2.15) and (2.16) on the right hand side of (2.13) one can now express / in terms of T>(f) 



and its first two moments as 
1 



/ 



X + <Po 



j-v(f) + 



X 



V(f)dfi - 



V{f)fidfx 



(2.17) 



Now a Chapman-Enskog-like expansion is performed in |6J. Therefore, a small parameter e is 
introduced in order to scale the emissivity and opacity terms 



J 



£ 



X 



X 
e 



which are considered to be large. Thus the expansion is performed for small mean free paths A = eA 
with A = (x + 4>o — (f>i)~ l ■ Inserting this scaling in (2.17) leads to 

/ 



X + <Po 



e X 



V(f)df J ,\-3 f xc/> l \- / V(f)fxdfi 



-i 



(2.18) 
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If we expand / = /o + sfi + £ 2 /2 + • • • and collect the terms of the same power of e in ( 2.18 ), using 
the linearity (and formally the continuity) of P, we obtain the zeroth order term 



fo 



x + 



1 + 



X 



I 

X 



and the first order term 



Bfl 



-1 



X + <f>o 



V{h) + ^ 



X 2 



P(/o)^ + 3^iA^ / P(/ )/id/i 



These expressions are now used to compute the approximation 



1 



the angularly integrated left hand side of (2.13), which by (2.14) is equal to j — x 



1 



fdfj, , the so- 



called particle number exchange rate with matter or total interaction rate. In order to simplify the 
computation, it is helpful to decompose the operator D additively into a part T> + that is symmetric 
with respect to p, and a part T>~ that is antisymmetric with respect to fx, i.e., we write 



V(f) = V + (f)+V-(f), 



with 



v + {f) 



cdt 



+ 



d In p 3v 

H 

cdt cr 



dp 



+ 



2 I din p 3v 
' cdt cr 



v 
cr 



CO 



doj 



and 



,5/ 
<9r 



^tt + -(1- A* 



2^ 
dp 



(2.19) 
(2.20) 



Calculating 



1 V(f +efr)dp = \ 
l ^ 



(2.21) 

we immediately see g /_x "D (fo)dfi = since V is antisymmetric in /i and /o does not depend on 
p. For the third summand we obtain 



V + {eh)dpi 



-i 



X + <Po 



p+(/ )+2?-(/o)+^ 



(2?+(/o)+2?-(/ ))d/i+3^iA- / (P+(/ )+P-(/o))^/i 



Since /q does not depend on /i, the term P + (/o)/x in the second inner integral of this expression 
is an antisymmetric polynomial in pi and thus vanishes after integration. By the same reasoning, 
the term D~(fo)p in the second inner integral is a symmetric polynomial in p so that angular 
integration gives an expression that no longer depends on p. Consequently, the last summand in 
the brackets is linear in p so that the application of the operator P + ( (•)) to this expression is 
an antisymmetric polynomial in p that vanishes after the outer integration. By the same reasoning, 
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the second summand in the brackets vanishes after the application of this operator and the outer 
integration. With these arguments and | f\ T>~(fo)d[i = as seen above we obtain the simplified 
equation 



1 * 



-1 



-i 



X + 



(2.22) 



in which the right hand side only contains terms that undergo the application of the operator T> + 
twice. In [6] these expressions are neglected in the calculation of s because they are considered 
"of higher order" than the "leading order term" ^ J^ 1 T> + (fo)dfj, whose contribution is already 
considered in (2.21). 

The last summand in (2.21) is given by 



V-{eh)dfi 



-i 



-i 



x + 



V + (f )+V-(f )+^l 

X 2 J- 



(P+(/o)+P-(/ ))^+3MaJ 
i z 



(P+(/o)+P-(/ ))Mi 



-1 



As for the third summand, the term T> + (fo)ii vanishes after integration in the second inner integral. 
Since /o is independent of fj,, the term T> + {fo), in the first inner integral is a symmetric polynomial 
in /j,. Therefore, since 

\ /_if(/o)d/i = the first inner integral does no longer depend on fi 
so that the application of the operator P~( (•)) to it is linear in jjL and vanishes after the 
outer integration. The same holds for the first term P + (/o) in the brackets which is a symmetric 
polynomial in fi so that the application of T>~ ( ^^ (•)) to it gives an antisymmetric polynomial in 
fj, that vanishes by the outer integration. Consequently, we obtain the equation 



V-{ehW=\ 
l * 



1 



-i 



x + 



p-(/o) + 3MA- / V-(f )fidfjL 



■•dfj,, 



that we can further simplify by observing 

2?"(/o) 

which leads to 
3/i0 



d_fo 
dr 



(2.23) 



1 f 1 1 

HA- J V (fo)ndn = SfifaX- 



1 ^ ~dr 



lA/x 



so that, with (2.12), we can conclude 

-l / 



1 V-(ehW = \ 
l z 



V 



i 



x + 



[l + (t>i\)V-(f 



\dfi 



dfo 
dr 



lAP-(/o), 



V 



-\Vr(J ))dn. 



Altogether we obtain the approximation 



\ J_ v+{hw - l - jT v- {xv- (/„)) dn , 



(2.24) 



of the total interaction rate s "in leading order". Note that there is no contribution to s in (2.21) 
that involves only one application of the operator T>~ . 
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To determine this approximation for s in concrete terms, we insert (2.19) and (2.20) in (2.24). 
First, we calculate 



V + (foW 



-I 1 


'dfo 


2 7-1 


cdt 


dfo + 


\ l ( 


cdt 


CO 



+ 



d In p 3v 

H 

cdt cr 



dlnp 
cdt 



cr J 



v 
cr 



.dfo 
duo 



+ 



d In p 3v 

H 

cdi cr 



cr 



.dfo 
du 



dfi 



dfo 1 dlnp ^dfp 
cdt 3 cdt <9cj 



Here, the second line is obtained by interchanging integration and differentiation, considering the 
notation (2.5) for fo in the first and third summand in the integral, using | /i 2 d/x = | for the 
latter and observing that the second summand vanishes since fo does not depend on fi. Recall 



that the same reasoning already led to the left hand side of the trapped particle equation (2.4). 



Therefore, here, the second term on the right hand side of (2.24) is the more interesting one that 



V- {\V-(fo))dn=l 
l z 



will lead to the definition of the diffusion source. Concretely, taking (2.23) into account, we see 



dfx 



1 

2 

ld_ 

3 dr 

ld_ 

3 dr 

r 2 dr 







dr 

dfo 
dr 

dfo 



A/i 



dr 



r 



+ 



1 



d_ 
d[i 

( M } ~dr~ M 



dfo 
dr 



-i 



dr 

2 A dfp 
3 dr 



3 r dr 



For the third and fourth line we used again \ /j, 2 d/j, = |, the isotropy of fo and A as well as the 
notation as in (2.5) for /q. The last line follows from the product rule. As a result, considering 



(2.7), we obtain a diffusion term induced by fo with A/3 as the (small!) diffusion parameter. 

Finally, for the derivation of E in (2.2) and ( |2.3| ) we consider / = /o + e/i in the Boltzmann 
equation (2.1) with the approximated collision integral as in (2.13) and set /* = fo and f s = ef\. 



Then, a first order approximation in e of the angularly integrated (2.1) or (2.13) gives 

,\df 



df Id hip df 
cdt 3 cdt doj 



^d_ 



3 dr 



) =j-x[f t + \f_J s d^j , (2.25) 



"in leading order", i.e., neglecting the terms in (2.22). On the right hand side we use again (2.5 



as well as (1.5). Note that f s does not need to be isotropic. Now, if we compare equation (2.25 



with the trapped particle equation (2.4) we obtain 

S diff := j— 1 <' 



+ 



X 



rdp. 



(2.26) 



-i 



3 9r J ' 2 

as a suitable definition for S in this limit. Since here, the first term is a diffusion term, equation 



(2.25) can be regarded as an approximation of the Boltzmann equation in the diffusion limit. 
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Therefore, we call £ a diffusion source. Observe from the above calculations that the diffusion 
term on the left hand side of (2.25) stems from the contribution of h J^_ 1 T>{f s )d^ to the total 
interaction rate "in leading order". We announce that in a forthcoming paper [I], the diffusion 
limit will be derived "in leading order" by a Chapman-Enskog expansion and, even without the 
"leading order" approximation, by a Hilbert expansion of the Boltzmann equation. 



The three coupled equations (2.4), (2.6) and (2.26) arise from taking the angular mean of (2.2) 
and (2.3) and (2.13), i.e., they are no longer dependent on /i. However, in spite of (1.5), the 



influence of the collision integral is contained "in leading order" in £ ^w,/*, \ Jj^ f s dfj,,j,J'j by 

its dependency on the mean free path A = (x + <fto — 

If the mean free path is too big, the diffusion limit ( |2.25 ) is no longer a good approximation of 
the Boltzmann equation (1.1) so that the diffusion term in (2.26) may become too big. Concretely, 
if we neglect the dynamics of the background matter, i.e., the second term on the left hand side in 



(2.4), then / tends exponentially towards the stationary state solution 



f 



s 



X 



(2.27) 



compare (4.2). Since > is an a priori condition for a distribution function, we obtain the 



necessary condition £ < j for the coupling term. In [6, p. 1177] this condition is obtained by physical 
arguments for large mean free paths where the streaming particle component f s dominates over the 
trapped particle component /*. In [1] the free streaming limit £ = j will be derived by a Hilbert 
expansion of the Boltzmann equation. 

Conversely, if the mean free path A tends to 0, i.e., e — > 0, then /* dominates over f s = ef\. In 
the limit e = we obtain £ = in (2.26). By (2.27) and the definition of x = j + Xi we have at 
least the requirement £ > — % since the distribution function always satisfies < 1. In [6, p. 1177] 
it is argued physically that should be at most jjx_i which was the zeroth order approximation 
of / above, because that function represents the distribution for thermal equilibrium. With this 
assumption one gets £ > from (2.27). In pQ the reaction limit £ = will also be derived by a 
Hilbert expansion of the Boltzmann equation. 

With these considerations regarding the free streaming and the reaction limits, the diffusion 



source in (2.26) is limited from above by j and below by in order to account for regimes where 



the diffusion limit is not valid. Altogether, with 



£ := 



mm < max 



1 



d^r 2 \df 

2 dr \ 3 dr 



f s dfi, 



-i 



(2.28) 



as the coupling term, the equations (2.4) and (2.6) with (2.5) and (2.8) give the Isotropic Diffusion 



Source Approximation (IDSA) of the Boltzmann equation as introduced in [B]. 



3 Reduced Boltzmann model problem 

In this section, we reduce the Boltzmann equation ( |1.1[ ) to a simpler one that will later serve as a 
model to perform first numerical tests. 



3.1 Assumptions 



For our reduction we decouple the Boltzmann equation (1.1) from the state of the matter that 



contributes to angular aberration ^ , Doppler shift emissivity j , opacity x an d scattering 



kernel R. The first assumption that we make is that the matter is frozen, so that its state does 
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not depend on time, which implies = = 0. We also assume that the emissivity and the 
opacity only depend on r, i.e, j(r) and x( r )- The third hypothesis is that there are no collisions 



and, therefore, R = and C(f) = 0. Under these hypotheses the Boltzmann equation (1.1) reduces 
to 

Qj. + V Q- r + ) =j(r)-x(r)f(t,r,n,u). (3.1) 

For simplicity, we rescale the time t' = ct and no longer mention the dependency on u from now 



the prime in the time scaling again, the reduced equation that we want to solve is 



on since equation (3.1) is monochromatic, i.e., oj does only appear as a parameter here. Dropping 

(3-2) 



df(t,r,n) df(t,r,fi) 1 2 df(t,r,n) 

Qj. + » + -(! - V ) ^ =j{r)~ X(r)f{t, r, //) 



3.2 Conservative formulation of transport equation 



In this subsection, we will show that equation (3.2) can be written in conservative form as 

df 



at 



+ v 



(3-3) 



where the divergence operator is given by V • 



Lemma 3.1 With the divergence operator just mentioned, we have V 



0. 



Proof 



By direct computation we calculate V- (iq ^ ^2^j = ?S r V + f ^(l-^ 2 ) = ^ - ^ = 0. 



Remark 3.2 The divergence operator we use can be derived from the 6D Cartesian divergence 
operator in phase space if we assume spherical symmetry and constant velocity. The radial term 
in the divergence operator is the usual radial term in spherical symmetry and the additional term 
represents the term for the velocity. 

3.3 Reduced form of the IDSA equations 

In order to compare the reduced Boltzmann equation with the equations of IDSA, we need to apply 



the same assumptions as above. The reduced system of equations corresponding to (2.4), (2.6) and 



(2.28) turns out to be the following. 



cdt 
2 rl 



1 d / r 

r 2 dr V 2 j_ 1 

1 d 



j - xf 
flidn 



X 



(3.4) 
(3-5) 



mm < max 



, 3 



The only changes compared with (2.4), (2.6) and (2.28) are that the Doppler shift term of the 



trapped particle equation is absent and that the definition of the mean free path now is A := \ . 
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4 Numerical treatment 



4.1 Grid 

The spherically symmetric computational domain for the IDSA, ^^> SA = [0, R] x [0,T], and the 
one for the Boltzmann model, fig ol = [0,R] x [—1,1] x [0,T], is represented by points (rj,i n ) and 
(ri,fij,t n ), respectively. We fix i max , j max and n max and define the grid by 

iR 2j n nT 

n+x/2 = -. , Hj+1/2 = -i + -. , t = - . 

'max Jmax "rnax 

We also define the cells dj by 

Cij = [r i - 1 /2,r i+l/2 \ x \p,j_y 2 , Vj+1/2] ■ 

4.2 Reduced model 
4.2.1 Time splitting 



In order to solve (3.3), we perform an order one time splitting by denoting 



t\M)---:i- \! mid t',(f):-Y-[ ( i (1 ^ 2) )/) ■ 



and writing equation (3.3), which is an autonomous ODE with respect to /, in the form 

^t = F(f):=F 1 (f) + F 2 (f). (4.1) 

Denoting the flow maps of the vector fields F, F\ and F 2 by $t,F, &t,Fi an d &t,F 2 i we approximate 
®t,F by ^tF '■= &tF2 $tFi- This is known as a Lie-Trotter splitting and gives an approximation 



of order 1 to the solution of (4.1), see [U p. 42]. 



4.2.2 Analytical treatment of the reaction term 

The flow map &t,Fi corresponds to the ODE 

at 3 XJ ' 



that has the analytical solution 



/(t) = /(0)e-*x + (l_ e -*x ) 4 ) (4 .2) 

X 



for every r 6 [0, R], fx £ [—1, 1] and uj £ [0,E]. It describes the exponential change from /(0) to 
the distribution function j/x which is known to be a thermal equilibrium function and, therefore, 
a Fermi-Dirac distribution 

3 1 



~ oj — 7 

X e~k¥~ + 1 



because of the fermionic nature of the neutrinos pp. 1177/9]. In the Fermi-Dirac distribution, 
7 is the chemical potential, k is Boltzmann's constant and 1? is the temperature of the matter. 
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-0.5 




/' 



0.5 



Figure 1: Flow map of the flux of neutrinos due to transport and without reactions. We can see 
that incoming neutrinos (// < 0) turn into outgoing neutrinos (/i > 0) and, therefore, eventually 
leave the computational domain again. 



4.2.3 Treatment of the transport term 

Figure [l] shows the behaviour of the flow map &t,F 2 m the domain (r, /x) G Q = [0, 10] x [—1, 1]. 

Remark 4.1 Figure [7] shows that the flow becomes large as the radius goes to zero. As a conse- 
quence, the CFL condition will be critical in the region where r is small if the overall flux F{f) is 
dominated by i^C/)- On the other hand, if in this region the overall flux is dominated by F\{f), 
then the exponential decrease towards equilibrium removes the dependency on /j, for small r as f 
becomes isotropic. 

We now derive a finite volume scheme for the transport term. The divergence theorem applied 
in our coordinates is given by 



V ( ^ J r 2 drdfi 



n T dn (r 2 drdn) , 



with TgQ^drdfi) denoting the trace measure. 

Integrating ( |3.3| ) on a cell Cij and defining := t^ttt f c f(t n , r, [i)r 2 drd[i with a forward Euler 
scheme, we obtain the finite volume scheme 



rn+1 rr, 

Jij — Jv 



At 



Mj + 1/2 f r i+l/2 



1 - ^+i/2)fifdr 



i-1/2 



+ 



1*3-1/2 



Mj + 1/2 



M/i+lj T i+l/2 



d\i 



n-1/2 

r i+l/2 



(l-M?-i/ 2 )/y-i^ 



for fj, < and 



f 



n+l 



Mj-l/2 ■''"i-1/2 



Mj-1/2 



n-1/2 



(l-M?-i/ 2 )/y-i^ 



n+1/2 



for /i > 0. This scheme uses a first order upwind approximation of the flux. 
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4.3 Implicit-explicit finite difference discretization of the IDSA 

In this section, we outline the discretization of the IDSA as introduced in {&, App. B] and adapt it 
to our case. 



4.3.1 Computation of the mean of the streaming particle 



The time integration of IDSA is done by taking equation (3.5), i.e., 

-1 



r 2 dr 



r 2 \ J f'fjdp 



1 



X: 



fdVL, 



1 



and integrating it over space 

1 

2 



r2 Jo 



fdfi 



r 2 dr . 



(4.3) 



There is no explicit time dependence in this equation as the streaming component is in the station- 



ary state limit, see Subsection 2.2 Hence, time integration of this component reduces to updating 
the angular mean of f s . Taking the zeroth and first moments of f s as variables, we discretize 



equation (4.3) explicitly by 

1 \ n+l 



ffidfi 



-1 



i+i 



E 



Xj 



fdn 



-1 



r ;( r ;+§ 



This is a mid-point rule formula. Once we have the flux 5/1 f s fid^L, we compute the updated 
I I-i Z^/" by discretizing (2.8) as 

2 



1 \ n + l 
J s d^ 



1 + 1 1 



max(rj,_R^) 



2_ 



max 



1 \ n+l 

f'lrin) A) 



As in P, p. 1188], we use the factor 



to convert the flux from the inner zone edge i— \ , where 



it has been computed, to the zone center i, where it is used in the computation of the diffusion 
source. The maximum operator forces the flux not to be directed against the density gradient. 
In [6, p. 1188], it is claimed that this increases the stability of the scheme. We did not test this 
statement here. The value of the neutrino scattering sphere radius i?" is considered as given in our 
model. 



4.3.2 Computation of the mean of the trapped particles 

In order to compute ^J.^dfi, we use equation (2.5) to simplify the notation. We start by 
discretizing equation (3.4), implicitly for stability reasons, by 



t,n+l 



f 

•1 1 



tji 



cAt 



■n+ 1 ~n+l f t,n+l v n+l 



which, by eliminating fj ,n+1 on the right hand side, can be rewritten as 



t,n+l 



A,n -n+l 
J i Ji 



cAt 



_ &n+l ft,n 
Xj J j 

l + x" +1 cAi 



n+l 



(4.4) 
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For ease of notation, we now define a := ^ ( — 7-2 ^t^t ) as the diffusion part of the diffusion 



3 dr 



source X. We discretize a by two centered finite differences as 



1 



2 \ /i+i /i „2 x . /' A'-l 
'„■ , 1 Aj , 1 - - I . i A, 



rj-ri-1 



3r f (ri+i-r^i) 
We write this expression more compactly as 

a l = -Wt + i-ft) + Ufi-fti), 
with the two quantities £j, £j defined by 

1 r i+i A *+i 1 



(4.5) 



6 



i) r i+ i 



Ci 



) - rj_i 



Consequently, the as yet unlimited written as is given by 

1 



with 



a: 



n 



n+l 



n+l 



< +1 + x" +1 



1 N, "+ 1 



er(A B i - /n + cm 



t,n\ 



/•t,n+lN 
/i-l J 



(4.6) 



(4.7) 



which is a semi-implicit discretization. We first perform all the computations without the limiters 
and apply the limiting only at the end of the computations. In this spherically symmetric example 
the diffusive fluxes propagate almost exclusively outwards, hence we choose to discretize explicitly 
the inward flux and implicitly the outward flux in (4.7). 



Inserting the equation (4.7) into equation (4.6) gives 

rin+l _ c n/ ft,n ft,n\ , j-ni ft,n+l ft,n+l\ . ~n+l 

— ~?i Ui+l - h ) + <>i Ui - Ji-l ) + X% 



1 s "+ 1 

f'diM 



(4. 



We compute the solution cell by cell from r = to r = R. Therefore, we can assume that all the 
i—1 indexed quantities are known in equation (4.8). The only term on the right hand side of (4.8) 



that is still unknown is /*' n+1 . To eliminate it we introduce the vanishing term — Cf/i'" + CP/j 



t.n 



in equation (4.8) that leads to 

v>n+l tnr ft,n ft,n\ . /-nt ft, 

^i — _ Si Ui+i - h ) + ^i Ui 



t,n+l\ , ~n+l 

T Xi 



1 N n + l 

f s dfl 



and eliminate the second term using equation (4.4) with the unlimited diffusion source We 
obtain 

-n+l ft,n _ 



■n+l 



'm+1 



QcAt 



Xj fj 
l + x7 +1 cAt 



Ki Ui+i-Ji J+C» Ui -Ji-i )+Xi 



1 N n + l 



Solving this equation for gives 



i + (cr +xD cAi 

+ (i+xr +1 cAt) 



CcAt(j 



n+l 



x? +1 /t n ) 



€i (/i+i - 7i ) + Ci Ui ~ /i-i ) 

n+l • 



rt,ra\ 



pi,n+l\ 
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We can now apply the limiter 



= min I max 



}■ 



The updated diffusion source can now be used in the computation of the updated f { ' n+ using 



equation (4.4) 



5 Numerical results 

In this section, we study the discrete reduced model and the discretized IDSA numerically. First, we 
perform a numerical convergence study in order to validate the discretizations and obtain concrete 
truncation error estimates for the study of the modeling error. Next, we illustrate the different 
behaviours of the two discretizations for a concrete example, for which we also quantify the modeling 
error of the IDSA. 

5.1 Study of the discretization error 

We study the discretization error in two different ways for the two models. In the case of the 
reduced Boltzmann problem, we pick a solution, define the corresponding j and x an d compare 
the computed solution with the exact one. In the IDSA case, however, this strategy does not 
work because we cannot find a simple solution of the coupled system. Therefore, we compare 
the solutions on different refinement levels with a reference solution on the most refined grid. 
Concretely, we choose 10 refinements of the coarse grid for the reference solution and, in order to 
compute the discretization order, compare the solutions obtained by 0, . . . , 8 refinements with it. 
The parameters used in this study correspond to those used in the numerical example shown later 
in Subsection 15.2.11 

To study the discretization error on the domain (r, //) 6 Q = [0, 10] x [—1, 1], we start from a 
coarse grid with N r = 5 discretization points in the r-direction, Nn = 3 discretization points in the 
/i-direction and with Nt = 1 time steps. We choose a time interval [0, T] where T = 0.0005 for the 
Boltzmann simulation and T = 0.1 for the IDSA one. 

Since we expect both discretizations to be of order one both in time and space, at each step we 
refine the grid by dividing the step sizes At, Ar and A/x by two. Therefore, these step sizes are 
proportional. The error is computed with the solutions obtained at the last time step. 

Concretely, in order to study the discretization error of the reduced Boltzmann model, we choose 
the analytical solution 

/ cxact :=^^[(r-10) 2 + (/,+ l) 2 + l] , 

which represents a paraboloid in the domain (t,r,(i) 6 [0, 1] x fi. The normalization constant 1050 
has been chosen to insure / cxact (r, /i, t) < 1. 

With the choice of the values for the emissivity j cxact and the opacity x oxact as 

• " "' :) {/4(r -10) 2 + (/,+ 1) 2 + 1]+ 2(1-^)^+1)} , 



u 1050 

and 



[(r - 10) 2 + + l) 2 + 1] - 2^(1 - t)(r - 10) 
Xcxact : ~ (l-t)[(r-10) 2 + Gu + l) 2 + l] 



the function / cxact is a solution of equation (3.3). 
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In this setting we now let T = 0.0005. By Remark |4.1| we need a very small value of At because 
we have At = T in the first iteration. It would be too costly to compute the order for T = 0.1 for 
this example because of the restrictive CFL condition. 

Figure [2] shows the discretization error in the infinity norm on Vt = [0, 10] x [—1, 1] for the 
reduced Boltzmann model in panel (a) and on [0, 10] for the IDSA model in panel (b). In both 
cases, the error is of order one as one would expect for simple linear cases with the discretizations 
used. 



Boltzmann Discretization error 








-2 -10 1 

10 10 10 10 



h 

(a) Reduced Boltzmann model with T = 
0.0005: The dependency of the error in the 
infinity norm on Q. on the spatial mesh size 
h = Ar is displayed. The curve with round 
markers represents the relative error, the curve 
with the plus markers displays the absolute er- 
ror. The straight line has slope one. 



IDSA Discretization error 








-3 -2 - 1 1 ^ 1 

10 10 10 10 10 



h 

(b) IDSA model with T = 0.1: The depen- 
dency of the error in the infinity norm on [0, 10] 
on the spatial mesh size h = Ar is displayed. 
In this case the absolute and the relative error 
are equal because the infinity norm of / is I . 
The straight line has slope one. 



Figure 2: Discretization error results for the reduced Boltzmann model and IDSA. 



5.2 Comparison of the reduced model and the IDSA 

In this section, we present a comparison of the reduced Boltzmann model and the IDSA by a 
numerical example and provide a study of the modeling error with the help of this example. 

5.2.1 Numerical example 

The numerical example we study is characterized by the following parameters. We use a domain 
(r, fx) 6 Bol = [0,R] x [—1, 1] with R = 9 and T = 20. The domain of computation r G n iDSA = 
[0, R] for the IDSA is different because we do not have the dependency on \x. The grid used here 
has N r = 100 discretization points in the r-direction and Nn = 20 discretization points in the ji- 
direction. The time step is At = 0.1. This time is small enough with regard to the CFL condition, 
because we force the region where the CFL condition is critical to be dominated by reactions, see 
Remark 14.11 

We choose the space interval to be [0, 9] and divide the computational domain into three parts 
in order to see the behaviour of the different (reaction, diffusion, free streaming) regimes of the 
Boltzmann equation that are representable in the IDSA. In our example, the occurrence of the 
different regimes depends on the values of the emissivity j(r) and of the opacity x(r) that are 
linked by a Fermi-Dirac distribution. We define the emissivity to be j(r) = 1 on [0,3], j(r) = 10 -3 
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on [6,9], linear on [3,6] and continuous on [0,9]. This choice leads to reaction dominance on 
[0,3] and dominance of transport on [6,9]. In the rest of the computational domain we have the 
transition between the two regimes. We fix the Fermi-Dirac distribution to be (exp(— 20) + 1) , 
which implies x = j (exp(— 20) + 1) for the opacity. Finally, we set the neutrino scattering sphere 
radius R u = 4.5. 

We use homogeneous Neumann boundary conditions at r = if \i < and at r = 9 if \i > 
and homogeneous Dirichlet boundary conditions at r = if \i > and at r = 9 if \i < 0. These 
boundary conditions are set to match the flow of transport. Dirichlet conditions correspond to 
incoming flow and Neumann conditions correspond to outgoing flow, see Figure [TJ The Neumann 
conditions are discretized by creating ghost cells just outside the domain that take the same value 
as the previous ones. For the IDS A we use the fact that we impose different regimes in different 
parts of the domain. We therefore use homogeneous Neumann boundary conditions at r = for /* 
and at r = 9 for f s and Dirichlet boundary conditions for the other two : f s (0) = and /*(9) = 0. 



Remark 5.1 We do not need any boundary conditions in the ^-direction because the flux vanishes 



along these boundaries, see Subsection J^.2.3 



As initial condition we use f(r, /U, 0) = \ for the reduced Boltzmann model. Motivated by the 
choice of j(r) we choose /* to be /*(r, 0) = \ on [0,3], /*(0) = on [6,9], linear on [3,6] and 
continuous on [0, 9] as well as f s (r, 0) = \ — /*(r, 0). 



Remark 5.2 A constant initial condition is not physical. In a spherical domain, we expect a 
decrease of the distribution function f with respect to the radius r. However, this is not a big issue 
since the reaction part will force any solution to converge to the correct one and the transport part 
will eventually remove the rest of the unphysical initial condition. 

The results of the simulation are shown in Figure [3] which displays four relevant snapshots 
of the evolution of the distribution function (panel (a)) and of the angular mean of it compared 
to the IDS A solution (panel (b)). Panel (a) of Figure [3] shows the distribution function at four 
chosen times. In the subdomain [0,3] x [—1,1], the evolution is driven by the reactions and, 
therefore, the only effect that we see is an isotropic growth of the function towards the equilibrium 
feq = j/x = (exp(— 20) + 1) _1 w 1. In the subdomain [6,9] x [—1, 1], the evolution is dominated 
by transport. It is therefore not isotropic and the neutrinos are moved along their trajectory lines 
as shown in Figure [TJ As explained in Remark 5.2 the unphysical initial condition is eventually 



removed. As the flow of transport does not have any effect on an isotropic zone because it is 
conservative, we start to see its effect on the boundary at r = 9 and for \i < 0. The fact that \x 
is negative reflects incoming neutrinos. As we set the distribution function to be zero outside the 
domain, the isotropy reduces the distribution function. The transport propagates in the subdomain 
as time evolves and removes the unphysical initial condition after a finite time. Another anisotropy 
is propagating in this domain, arising from the center reaction dominated zone. This zone acts 
as a source of neutrinos that propagate outwards, into the [i > subdomain. We start to see 
this effect in the t = 2 snapshot, and it continues to grow in the other two. At the end of this 
simulation, the unphysical initial condition has been completely removed, the reaction dominated 
region is equilibrated and produces neutrinos that propagate mainly outwards, driven by the flow 
of transport. The intermediate region just shows how the reaction dominated region is coupled 
with the transport dominated region. 

Panel (b) of Figure [3] shows the comparison between the two models at the same four times 
as panel (a). In order to discuss the comparison, we divide the domain into the 3 subdomains as 
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(a) Boltzmann distribution function in the phase space 
(r,/i) 6 [0,9] x [-1,1]. 



(b) Comparison of the averaged distribution functions. 
The solid blue line represents the reduced Boltzmann 
solution. The blue dashed line represents the IDSA 
solution, which is the sum of the red dashed line repre- 
senting the trapped component and the green dashed 
line representing the free streaming component. 



Figure 3: This figure shows 4 snapshots of the evolution of the Boltzmann distribution function / 
in the phase space (panel (a)) and the same for its angular average compared to the IDSA solution 
(panel (b)) at times t = 0.7, 2, 7.5 and 20. 

before, concretely into f2 rcac = [0,3], the subdomain dominated by reactions, transp = [6,9], the 
subdomain dominated by transport and the intermediate subdomain O tnt = (3,6). 

In O roac we see that the convergence of the IDSA solution / IDSA to the equilibrated one is much 
faster than the convergence of solution / Bol for the Boltzmann model. This is explained by the 
fact that the trapped component of /idsa is strongly coupled with the free streaming component 
which evolves infinitely fast since we use a stationary state approximation. The coupling between 
the two components of /idsa therefore explains the faster convergence rate of it to equilibrium. As 
seen in the discussion of panel (a), in this domain, the distribution function / Bol is isotropic and, 
therefore, its angular mean evolves in the same way. We notice that the IDSA is underestimating 
the Boltzmann solution. 

In f^transp we see that the angular mean of / Bo i describes the reduction of the unphysical initial 
condition as explained in the discussion of panel (a). In the two first graphs at times t = 0.7 
and t = 2, we see some regions that are still representing the isotropic initial condition. In the 
two other graphs, we do not notice its presence anymore even if there is still a component of it 
at t = 7.5 as shown in panel (a). For the IDSA we see, as expected, that this region is described 
by the streaming component of / IDSA , but it shows an underestimation of the reference Boltzmann 
solution. We think that this should be explained by the way of coupling the two regions rj reac and 

^transp • 

In the intermediate subdomain tt int , we see the evolution of the transition between the two 
other subdomains, in particular, the transition between the trapped and the streaming components 
of /idsa- As expected, we see the growth of the streaming part in the intermediate domain and a 
decrease in the transport domain. There are two aspects that we want to point out here. First, 
we see a strange behaviour of the streaming component around r = 4.5. This behaviour is a 
consequence of the neutrino scattering sphere radius that has been set precisely to this value of 
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r. Second, we observe that the streaming component does not vanish in the reaction dominated 
regime as expected. There is a noticeable component until r ~ 1 and, as a consequence of the 
coupling, this reduces a little the trapped component which, then, underestimates the equilibrium 
value. 



5.2.2 Study of the modeling error in IDSA 



In this subsection and on the basis of the numerical example of the previous subsection, we study 
the modeling error of the isotropic diffusion source approximation compared to the Boltzmann 
model. In order to have a reliable measure of this error, we need to control the discretization error. 



Using the results of Section 5.1 we choose the grid parameters such that the relative discretization 
error is less than one percent. The parameters are iV r = 257, Np = 129, and At = 1.5625 • 1CP 3 . 
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(a) IDSA modeling error, computed as the relative error (b) Localization of the error at time t=100, computed by 
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Figure 4: Modeling error of the IDSA. 

The results shown in Figure [4] quantify the error of the IDSA with respect to the Boltzmann 
model. As we fix the discretization error to be smaller than one percent, we know that it contributes 
less than 0.02 in the graph of panel (a). We can therefore conclude that the relative modeling error 
is 20% ± 2% which one might find quite big. A relevant question to ask is if the error is uniform 
or accumulates in some regions. To answer this question, we compute in each spatial point the 
absolute difference between the two solutions. The result is shown in panel (b) of Figure |4j As 
expected, at the end of the simulation, the error is mainly located in the region of transition (3, 6) 
between the reaction and free streaming regimes. 

Summarizing, this preliminary study shows that the IDSA is qualitatively reasonable, but it 
also shows that the coupling between the two regimes exhibits considerable errors and is the main 
source of error of the IDSA. 



References 

[1] H. Berninger, E. Frenod, M. J. Gander, M. Liebendorfer, J. Michaud, and N. Vasset. Derivation 
of the IDSA for Supernova Neutrino Transport by Asymptotic Expansions. In preparation, 
2012. 



21 



[2] S. W. Bruenn. Stellar Core Collapse: Numerical Model and Infall Epoch. ApJS, 58:771-841, 
1985. 

[3] R. Buras, H.-T. Janka, M. T. Keil, G. G. Raffelt, and M. Rampp. Electron Neutrino Pair 
Annihilation: A New Source for Muon and Tau Neutrinos in Supernovae. ApJ, 587:320-326, 
2003. 

[4] E. Hairer, C. Lubich, and G. Wanner. Geometric Numerical Integration. Springer, 2002. 

[5] M. Liebendorfer, O. E. B. Messer, A. Mezzacappa, S.W. Bruenn, C.Y. Cardall, and F.-K. 
Thielemann. A Finite Difference Representation of Neutrino Radiation Hydrodynamics in 
Spherically Symmetric General Relativistic Spacetime. ApJS, 150:263-316, 2004. 

[6] M. Liebendorfer, S. C. Whitehouse, and T. Fischer. The Isotropic Diffusion Source Approxi- 
mation for Supernova Neutrino Transport. ApJ, 698:1174-1190, 2009. 

[7] A. Mezzacappa and S. W. Bruenn. Type II Supernovae and Boltzmann Neutrino Transport: 
The Infall Phase. ApJ, 405(2):637-668, 1993. 

[8] D. Mihalas and B. Weibel-Mihalas. Foundation of Radiation Hydrodynamics. Oxford Univer- 
sity Press, 1984. 

[9] M. Rampp and H.-T. Janka. Radiation hydrodynamics with neutrinos - Variable Eddington 
factor method for core-collapse supernova simulations. Astron. Astrophys., 396:361-392, 2002. 

[10] E. T. Whittaker and G. N. Watson. A Course of Modern Analysis. Cambridge University 
Press, 4th edition, 1980. 



22 



